\documentclass[12pt]{article}

\usepackage{HWStyle}


\newcommand{\hmwkClass}{Programming Task \#1}
\newcommand{\hmwkAuthorName}{李惠腾}
\newcommand{\hmwkTitle}{Design Document}

 \setlist[3]{noitemsep}

\begin{document}

\maketitle

\section{问题介绍}
数值解课程讲义的第7章介绍了拉普拉斯算子的二阶离散格式.
其中,对于Dirichlet条件的不规则边界问题,我们给出了一套
在正则内部点有二阶精度,在近边内部点有一阶精度的离散策略,并利用
离散极值原理(DMP)证明了它在max-norm范数下的全局二阶收敛性.

这个启发性的算法框架带给了我们两个方面的新问题. 
第一,在理论上,我们这套混合精度算法是否能推广到其他边界类型的椭圆算子问题上去？对于Neumann边界或是混合边界类型的问题,我们是否也能选择一套在正则内部点有高精度,
在近边内部点稍低精度的算法,使之具有给定精度的全局收敛性？
第二,如果理论上证明可行,那如何设计一个尽可能简便且二阶收敛的混合精度框架？
本文将在附录A中简要回答第一个问题, 而在正文用一套具体的算法设计框架来回答第二个问题. 

不失一般性,考虑计算域$\Omega$是一个内部挖去圆盘$D$的正方形($D$的半径可为0). 在$\Omega$内,它被$-\Delta u = f$控制. 而在方形的四条边与$\partial D$上,边界条件
被Dirichlet(`D')或Neumann(`N')中的其一控制,有31种有限的组合(不考虑纯Neumann边界条件). 

\section{算法原理}
对于方形边界和圆形边界,我们均采用ghost cell离散策略. 

这给编程带来的直接好处就是不用再受制于近边内部点的复杂求解格式. 
我们可以在近边内部点上直接使用标准五点差分格式, 而把上述过程所需、但在计算域之外的网格点
作为一个待求解量也加入方程组, 并利用合适的边界条件来填充这些\textbf{ghost cell stencil (GCS)}. 
这种边界与内部解耦的策略能带来的另一好处是可以自由控制在边界附近离散的精度,从而更轻松地实现混合精度框架. 

本次程序所需要的共有四种GCS.  

\begin{enumerate}
    \item 方形边界 Dirichlet条件
          
          对于边界上的点$P$, \[U_P = u_P.\] 显然误差为0.
          
    \item 不规则边界 Dirichlet条件
          
          \begin{figure*}[!htbp]
              \centering
              \includegraphics[width=0.45\textwidth]{cirDiri.jpg}
              \caption{The points used for discretizing the Dirichlet boundary condition.}
          \end{figure*}
          
          以讲义中的图1为例,对于计算域外的点$N$, 我们用线性插值或三点插值来填充GCS:
          \begin{itemize}
              \item \begin{align}
                        \alpha U_{N} + (1-\alpha)U_{P} = u_{B}.
                    \end{align}
                    线性插值$u_B$的局部截断误差\textbf{(LTE)}为$O(h^2)$.
              \item \begin{align}
                        \frac{(\alpha+1)\alpha}{2}U_{N} -(\alpha+1)(\alpha-1)U_{P}+\frac{(\alpha-1)\alpha}{2}U_S=u_B.        
                    \end{align}
                    三点插值$u_B$的LTE为$O(h^3)$.
          \end{itemize}
          
          求解器会同时实现这两种并将选择权作为一个预置选项. 
          
          首先, 方案(2)的本质就是讲义上提供的一阶精度不规则格式(7.86). 因为当我们在点$E$也用三点插值,
          再将它们与$P$处的五点标准差分格式结合在一起后, 消去$U_N$,$U_E$便是LTE为$O(h)$的(7.86).
          
          其次, 我们其实可以直接选用第一个公式. 虽然它结合上五点差分格式后的LTE会
          掉到$O(1)$, 但这时定理7.63中的$U_A$
          的系数为$-\frac{1}{\theta h^2}<-\frac{1}{h^2}$同样成立,
          且$T_2=K_2$时, 也
          能有$|E_P|\leq (F_1R^2+F_2)\max\{\frac{K_1}{4F_1},\frac{K_2}{2F_2}\}h^2$, 
          故方案(1)也是能使终解误差二阶收敛的. 
          
    \item 方形边界 Neumann条件
          
          对于边界上的点, 我们仍使用标准五点差分格式. 于是边界外围会生成一层ghost cell. 
          对于ghost cell $P$, 我们用在紧挨$P$的边界网格点$Q$上中心差分公式来填充GCS,
          \begin{align}
              \frac{1}{2h}U_P - \frac{1}{2h} U_R = g_Q.     
          \end{align}
          
          其中$R$是一个满足$Q$为$PR$中点的内部点, $g_Q$为外法向导数在$Q$的值. 
          用中心差分公式来逼近$\frac{\partial u}{\partial n}(Q)$时, LTE为
          $O(h^2)$.
          值得注意的是,当我们把这个GCS与对应的五点差分格式耦合时,得到的四点公式是只有一阶
          LTE的.
          
          
    \item 不规则边界 Neumann条件
          
          \begin{figure*}[!htbp]
              \centering
              \includegraphics[width=0.6\textwidth]{cirNeum.jpg}
              \caption{The points used for discretizing the Neumann boundary condition.}
          \end{figure*}
          
          以图2为例,对于计算域外的点$P$, 我们提供两种填充不同精度的GCS策略. 
          
          首先,自圆心向$P$引一条射线. 由于曲线是圆,故射线的反向就是$J$处外法方向. 
          通过泰勒展开,我们能得到外法向导数的一阶LTE近似公式,
          \[\frac{\partial u}{\partial n}(J) = \frac{u_P-u_{P_1}}{L}+O(h),\]
          其中$L:=PP_1=\frac{h}{\sin(\alpha)}$.  而对于$u_{P_1}$,我们使用$M_1,C_1$上的线性插值公式(我们把$M_1$叫作\texttt{faceMid}, $C_1$叫作\texttt{Corner}),
          \[
              u_{P_1}=(1-r)u_{M_1}+ru_{C_1}+O(h^2), 
          \]
          其中$r:=\frac{M_1P_1}{h}=\cot(\alpha)$,是$P_1$相对于\texttt{FaceMid}的截距比. 将两个式子合并起来就是$P$处的GCS,
          \begin{align}
              \frac{1}{L}U_P-\frac{r}{L}U_{C_1}-\frac{1-r}{L}U_{M_1} = g_J,
          \end{align}
          $g_J$为外法向导数在$J$的值. 此时LTE为$O(h)$. 
          
          这个一阶算法有两大好,一是确保了耦合方程下的DMP一定成立. 由于充分利用了曲线几何信息来确定GCS的选取,我们能保证$r\in[0,1]$,
          故不会给等价的耦合方程组带来非对角占优. 二是确定GCS相当简便,尤其是对于等距网格. 
          事实上,对于任意$\alpha$, $\alpha$为射线与x轴(逆时针为正向)的夹角 ($\alpha \in [-\pi,\pi]$),我们都有
          $$L=\frac{h}{\max(|\sin\alpha|,|\cos\alpha|)},\quad r=\min(|\tan\alpha|,|\cot\alpha|).$$
          另外在选取上,只要知道$\alpha/(\pi/4)$的值,便可唯一确定\texttt{FaceMid}与\texttt{Corner},实现很简单. 
          
          当然,它的坏就是当它耦合
          上标准五点差分模板后, LTE会只有$O(1)$. 与不规则Dirichlet边界不同, 这时$O(1)$是不够用的. 
          事实上, 若启用了它,在理论上我们只能证明求解框架是全局一阶收敛的. 并且数值实验也显示了这一现象.
          
          作为补救,我们其实也能把这个方法延拓到二阶,只要我们将$P_2$也引入外方向导数估计中,
          \[
              \frac{3L-L_1}{2L^2}U_P - \frac{2(L-L_1)}{L^2}U_{P_1}-\frac{2L_1-L}{2L^2}U_{P_2}=g_J,  
          \]
          其中$L_1:=PJ$.接着在$P_i$上用三点插值,剩余的具体细节与一阶方法类似. 值得注意的是, 这时我们无法确保DMP必然成立,因为$2L_1$不是总会大于$L$. 
\end{enumerate}

\section{程序设计}
\subsection{基本类型}

\begin{itemize}
    \item \texttt{\#include<Vec.h>, \#include<Tensor.h>, \#include<Box.h>}
    \item 
          \texttt{using Real = double;\\
              enum\{Dim = 2\};\\
              using iVec = Vec<int, Dim>;\\
              using rVec = Vec<Real, Dim>;\\
              template <class T>\\
              using vector = std::vector<T>;
          }
    \item 
          \texttt{enum directionStyle \{ \\
              Clockwise = -1,\\
              CounterClockwise = 1\}  
          }
\end{itemize}

\subsection{class OrientedCircle}

\begin{itemize}
    \item 表示平面上的有向圆(曲线定向的左侧为曲线内部).
    \item 数据成员
          \begin{verbatim}
protected:
    Real radius;
    rVec center;
    directionStyle direction;
    Real tolerance;
\end{verbatim}
    \item 函数成员
          \begin{enumerate}
              \item \texttt{OrientedCircle(const rVec \&aCenter, const Real aRadius,\\
                        const directionStyle aDirection = Clockwise, \\
                        const Real atol = 1e-15) \\
                        : center(aCenter), radius(aRadius), direction(aDirection),tolerance(atol)\{\}}
                    \textbf{Input:} 圆心坐标,半径,定向与不确定性因子.  \\
                    \textbf{Effect:} 初始化参数. 默认圆心在曲线外,不确定性控制因子为$1e-15$.
                    
              \item 基本操作(幅角与圆上点的互相转化,计算圆上给定点的外法向量,计算能覆盖圆的最小方形的四个顶点)
              \item \texttt{bool inEnclosure(const rVec \&pos) const}\\
                    \textbf{Input:} A point \(pos\) of type rVec. \\
                    \textbf{Output:} 返回true说明该点位于曲线围成区域的内部,返回false为外部或边界.
                    
              \item \texttt{const vector<rVec> intersection2CartesianLine(int d, Real pos\_d) const} \\
                    \textbf{Input:} 输入d方向分量的数值x\_d. \\
                    \textbf{Output:} 计算并返回$x_d=$pos\_d与圆的交点.
              \item \texttt{const rVec intersection2CartesianLine(int d, const rVec \&pos) const }
                    \textbf{Precondition:} $x_d=pos[d]$与圆的交点非空. \\ 
                    \textbf{Input:} 输入一个d方向与x点. \\ 
                    \textbf{Output:} 计算$x_d=pos[d]$与圆的交点,并返回其中距离pos点最近的交点.
          \end{enumerate}
          
\end{itemize}

\subsection{class RectDomain}

\begin{itemize}
    \item 表示二维平面上的矩形网格,一个被赋予了距离的\texttt{box}继承类.
    \item \texttt{template <int Dim>\\
              class RectDomain : public Box<Dim>}
    \item 数据成员\\
          \texttt{public:\\using Box<Dim>::lo;\\
              using Box<Dim>::hi;\\
              protected:    rVec dx;}
    \item 基本操作(返回网格间距,判断实坐标是否落在矩形中)
\end{itemize}
\subsection{class ExcludedCircle}

\begin{itemize}
    \item 表示内部嵌入了一个圆盘的矩形网格. 额外记录了圆盘与网格之间的互动信息.
    \item \texttt{template <>\\
              class ExcludedCircle<2> : public RectDomain<2>}
    \item 数据成员\\
          \texttt{public:\\
              const OrientedCircle embeddingCircle;\\
              const Real tol;\\
              Tensor<bool, Dim> inDomain;///<存储了网格上每个点是否在计算域中.true说明在内部,false说明在外部或边界.\\
              const Box<2> getMiniBox() const; ///<存储了能包含圆的最小网格盒.\\\\
              ///<代表两个不同方向上的截断比矩阵. \\
              ///<第d个Tensor的(i,j)坐标存储了从（i,j）沿e\_d正方向出发到intersectionPoint的相对网格距离比例.若0,则说明无交点. \\
              ///<这里假定每条mesh segment至多与圆只有一个交点. \\
              ///<于是我们可以把每个交点向下/左投影到最近的网格节点,并在这些网格点存储对应线段的截断比$\alpha$.\\
              ///<其中$\alpha$代表了节点（关于下/左）的相对比例. \\
              Tensor<Real, Dim> cutRadios[2];\\\\
              ///< 将截断比矩阵中的截断比转化为实坐标\\
              const rVec cutPointMap(iVec \&idx, int d) const;
          }
    \item \texttt{ExcludedCircle<2>::ExcludedCircle(const RectDomain<2> \&aRectDomain,\\
              const OrientedCircle \&aCircle, Real atol)\\
              : RectDomain<2>(aRectDomain), embeddingCircle\{aCircle\}, tol(atol);
          }\\
          这个类上唯一且最重要的操作是唯一指定的构造函数. 在构造过程中,我们初始化并先后确定\texttt{miniBox}, \texttt{inDomain}与\texttt{cutRadios}.
\end{itemize}

\subsection{class DirectPossionSolver}

\begin{itemize}
    \item 用于求解目标问题的求解器. 为计算域中的点和ghost cell 编号,组装系数矩阵. 输入右端项数据后,直接LU分解求解,误差分析与格式化输出.
    \item \texttt{template <>\\
              class DirectPossionSolver<2>;}
    \item 数据成员\\
          \texttt{protected:\\
          ExcludedCircle<2> emGrid;\\
          ///<记录边界类型. 只支持'D'或 'N', 其中前四个为方形边界类型,对应关系是$d$方向的$side$边对应了bcTypes[2*d+(side+1)/2], 其中 $d\in\{0,1\},side \in\{-1,1\}$.第五个是内部圆形边界类型. \\
          char bcTypes[5];\\
          ///<记录方程未知量中计算域里的所有点(包含方形边界上的点)的总数.\\
          int numUnknownsInDomain;\\
          ///<记录方程未知量中除去不规则边界上的ghost cell的总数. 即,计算域里的所有点,加上方形边界上的所有Neumann边界对应的ghost cell的总数.\\
          int numUnknownsExceptIrGhost;\\
          ///<记录方程全体未知量的总数.\\
          int numUnknownsAll;\\
          ///<索引映射. 输入格点整数坐标,输出作为方程未知量编号. 若不是方程未知量,返回-1.\\
          Tensor<int, 2> unknownMap;\\
          ///<索引逆映射. 输入作为方程未知量编号,输出格点整数坐标.\\
          vector<iVec> unknownMapInv;\\
          ///<系数矩阵\\
          Tensor<Real, 2> poissonMat;\\\\
          ///<控制不规则边界ghost cell stencil的精度.\\
          int ordIntpDiriIrBdr = 1;\\
          int ordIntpNeumIrBdr = 2;\\\\
          public:\\
          ///< 唯一的构造函数. 输入嵌入网格与边界类型. 效果为生成索引(逆)映射和系数矩阵.\\
          DirectPossionSolver(const EmGrid \&aGrid, const char *aBcTypes);\\\\
          ///< 求解方法1: 输入右端项信息(规则区域已经离散完成). 输出数值解.\\
          template <class TFunc>\\
          Tensor<Real, 1> computeSolution\\
          (const Tensor<Real, Dim - 1> *aRectBdrData,\\
          const Tensor<Real, Dim> \&aDomainData,\\
          const TFunc \&aIrBdrDataExpr) const;\\
          ///< 求解方法2: 输入右端项函数信息. 输出数值解.\\
          template <class VecExpr, class ScalerExpr>\\
          Tensor<Real, 1> computeSolution\\
          (const ScalerExpr \&negDDexpr,\\
          const VecExpr \&Dexpr,\\
          const ScalerExpr \&expr) const;\\
          }
    \item 
          \begin{verbatim}
DirectPossionSolver<2>::DirectPossionSolver(const EmGrid &aGrid,
                                            const char *aBcTypes)
: emGrid(aGrid)
{
strncpy(bcTypes, aBcTypes, 5);
const auto dx = emGrid.spacing();
// Step 1. Find 1-to-1 correspondence between mesh points and unknowns of matrix.
doFindUnknownMap();
// Step 2. Fill in the coefficients of matrix.
doFillMat();
}
\end{verbatim}
          我们将在构造函数中先后生成索引映射和系数矩阵. 特别地,对于\texttt{doFillMat},我们的处理步骤如下:
          \begin{enumerate}
              \item    Fill \texttt{possionMat} for the regular ghosted boundary.
                    For `D' we fill the boundary and for `N' we fill the ghost.
                    Also generate the box for all points of 5-point-standard-stencil.
                    
              \item Fill \texttt{possionMat} for the standard box, i.e. points that support 5-point-standard-stencil.
                    
              \item Fill \texttt{possionMat} for irregular boundary part.
          \end{enumerate}
    \item 对于\texttt{computeSolution}, 我们通过重构提供了两种输入方案. 一种是输入已离散好的规则区域及边界的Tensor以及不规则边界的
         函数信息, 第二种是直接输入所有函数信息. 
    \item  误差分析与格式化输出
          \begin{verbatim}
public:
/**
* @param  q       q-norm of grid function. max-norm when q = 0.
*/
template <class TFunc>
const Real computeErrNorm(const TFunc &exactSolution, const int q = 1) const;
/**
* @brief 将计算域网格上的数值解,真解与误差以方阵形式打印到指定文件中. 
* 空缺部分用0占位. 
* @param  _file  给定文件名.          
*/
template <class TFunc>
const void plotIn(const std::string &_file, const TFunc &exactSolution) const;  
\end{verbatim}
\end{itemize}

\subsection{class FuncFiller}

\begin{itemize}
    \item 这个类主要用来将给定的函数填充到离散网格点上.
    \item 数据成员\\
          \texttt{protected:\\
              RectDomain<Dim> rd;
          }
    \item 基本操作:填充计算域内值,方形边界值以及方形边界外法向导数.
          \begin{verbatim}
public:
template <class TFunc>
void fillDomainVal(Tensor<Real, Dim> &aData, const TFunc &negDDexpr) const;

template <class TFunc>
void fillRectBdrVal(Tensor<Real, Dim - 1> &aData, int D,
                   int side, const TFunc &expr) const;

template <class TFunc>
void fillRectBdrNormalGrad(Tensor<Real, Dim - 1> &aData, int D, 
                  int side, const TFunc &Dexpr) const;    
\end{verbatim}
\end{itemize}

\newpage
\section{附录A}
第七章中定理7.63的证明过程仅适用于例子7.59这种纯Dirichlet边界问题. 而对于这次项目所要处理的
混合边界问题, 这里以一种情况为例加以补充证明. 

\begin{figure*}[!htbp]
    \centering
    \includegraphics[width=0.6\textwidth]{DMPexp1.jpg}
    \caption{}
\end{figure*}

如图3, 对内部挖去圆盘的$[0,1]\times[0,1]$计算域, 考虑这样一种边界情况: 在四个方形边界都是Neumann条件, 在内部圆形边界是Dirichlet条件. 
接着我们
在计算域五类典型点处给出计算模板.
\begin{enumerate}
    \item 对于东侧点$P_1$, 
\begin{align}
    L_hU_{P_1}:=\frac{4U_{P_1}-2U_{W_1}-U_{N_1}-U_{S_1}}{h^2},
\end{align} 
对应的右端项为$f_{P_1}+\frac{2g_3(P_1)}{h}$. 这个等式引入的LTE为$O(h)$. 
\item 对于东北角的点$P_2$,
\begin{align}
    L_hU_{P_2}:=\frac{4U_{P_2}-2U_{W_2}-2U_{S_2}}{h^2},
\end{align}
对应的右端项为$f_{P_2}+\frac{2g_3(P_2)+2g_1(P_2)}{h}$. 这个等式引入的LTE为$O(h)$. 
\item 对于正则内部点$P_0$,
\begin{align}
    L_hU_{P_0}:=\frac{4U_{P_0}-U_{W_0}-U_{N_0}-U_{S_0}-U_{E_0}}{h^2},
\end{align}
对应的右端项为$f_{P_0}$. 这个等式引入的LTE为$O(h^2)$. 
\item 对于靠近圆的近边内部点$P_3,P_4$, 我们将之前的GCS方案(1)与标准五点差分模板耦合起来作为计算模板, 引入的LTE为$O(1)$.
\end{enumerate}
下面我们证明这套混合精度求解框架确实在max-norm下二阶收敛.
\begin{proof}
    首先定义试验函数
    \begin{align}
        \Phi(P)=\left\{\begin{array}{lc}(x-2)^2+(y-2)^2+1,&P\in X_{\partial\Omega}\\(x-2)^2+(y-2{)^2,}&P\in{\mathrm X}_\Omega\;\end{array}\right.
    \end{align}
    
    需要注意, 如例子7.54中所述, 方形Neumann边界上的所有网格点同样是 equation-discretization point, 即属于$X_\Omega$而非$X_{\partial\Omega}$.
    故$X_{\partial\Omega}$就是圆形Dirichlet边界. 接下来我们考虑$X_\Omega$中各个点处$L_h\Phi$的结果.
    
    \begin{enumerate}
        \item 对$P_1$, 
    \begin{align*}
        L_h\Phi_{P_1}&=\frac{4\Phi_{P_1}-2\Phi_{W_1}-\Phi_{N_1}-\Phi_{S_1}}{h^2}\\
        &=-4-\frac{4}{h}<-\frac{4}{h}.
    \end{align*} 
    \item 对$P_2$,
    \begin{align*}
        L_h\Phi_{P_2}&=\frac{4\Phi_{P_2}-2\Phi_{W_1}-2\Phi_{S_2}}{h^2}\\
        &=-4-\frac{8}{h}<-\frac{4}{h}.
    \end{align*}
    \item 对$P_0$,
    \begin{align*}
        L_hU_{P_0}=-4.
    \end{align*}
    \item 对$P_3,P_4$,
    \begin{align*}
        L_h\Phi_{P}\leq -4-\frac{1}{h^2}<-\frac{1}{h^2}.
    \end{align*}
    \end{enumerate}
    接着我们把$P_1,P_2$类型点的最大LTE记做$T_1=K_1h$, 把$P_0$类型点的最大LTE记做$T_0=K_0h^2$, 
    把$P_3,P_4$类型点的最大LTE记做$T_2=K_2$. 不难把定理7.61推广到$\mathbf{ X_{\Omega}=X_1\bigcup X_2 \bigcup X_3}$的情形, 
    且经验证$L_h$满足(DMP-1,2). 故我们有:
    \begin{align}
        \forall P \in \mathbf{X}, \; |E_P|&\leq \left(\max_{P\in\mathbf{X}_{\partial\Omega}}\Phi(Q)\right)\max\left\{\frac{K_0h^2}{4},\frac{K_1h}{4/h},\frac{K_2}{1/h^2}\right\}\\
        &<9\max\left\{\frac{K_0}{4},\frac{K_1}{4},K_2\right\}h^2=O(h^2).
    \end{align}
    故这套求解框架在max-norm下二阶收敛.
\end{proof}
可以直接观察到的一点是, 对于这个问题, 
在这些典型点任意混合精度求解框架的LTE阶数均不能低于证明中所使用的LTE阶数,
不然式子(9)将无法保有$O(h^2)$, 我们也自然无法用DMP来判定它二阶收敛. 
更广泛地说, 对于任意一个已经满足DMP条件的混合精度算法框架,
如果要让它在max-norm下$k$阶收敛即$|E_P|=O(h^k)$, 只要让:
\begin{enumerate}
    \item 在正则内部点的LTE至少$O(h^k)$.
    \item 在靠近Dirichlet条件边界的内部点的LTE至少$O(h^{k-2})$.
    \item 在靠近Neumann条件边界的内部点的LTE至少$O(h^{k-1})$.
\end{enumerate}

回到我们的目标问题. 对于其他类型的混合边值条件, 证明基本是类似的, 但除了带有圆形Neumann边界的问题. 
由于没有了规则边界时的对称性, 在估计外法向导数时候我们用不了中心差分格式, 
自然也无法享受到中心差分带来的额外一阶精度的LTE收益. 
这会直接导致的问题是使用GCS方案(4)时会把$|E_P|$降至$O(h)$. 
当然, 我们
自然可以使用更高阶的差分公式来估计边界法向导数, 但代价是会在GCS中使用更多的点. 而
通常这样的GCS耦合上标准五点差分模板后是没法满足(DMP-1,2)的. 对于本文中所使用的
估算外法向导数的第二种方法, 一方面, 我们无法确保DMP必然成立,因为$2L_1$不是总会大于$L$. 
另一方面， 数值实验又表明了它在绝大多数场合下的二阶收敛性. 这也表明, DMP的条件本质是一种基于对角占优的强稳定性条件.这种强稳定性对于
终解误差的收敛阶来说, 是严格且充分不必要的.   
\end{document}